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' Abstract 

' This paper presents the scientific body of knowledge behind the Human Biodynamics En- 

gine (HBE), a human motion simulator developed on the concept of Euclidean motion group 
<—{ • SE(3), with 270 active degrees of freedom, force-velocity-time muscular mechanics and two- 

h—i ' level neural control - formulated in the fashion of nonlinear humanoid robotics. The following 

QQ , aspects of the HBE development are described: geometrical, dynamical, control, physiological, 

' AI, behavioral and complexity, together with several simulation examples. 

I I , Index Terms: Human Biodynamics Engine, Euclidean SE(3)-group, Lagrangian/Hamiltonian 

' biodynamics. Lie-derivative control, muscular mechanics, fuzzy-topological coordination, 

^) \ biodynamical complexity, validation, application 

.2 ■ 

; 1 Introduction 

I i ' Both human biodynamics and humanoid robotics are devoted to studying human-like motion. 

They are both governed by Newtonian dynamical laws and reflex-like nonlinear controls [J [51 [3J 

^ ; aiSlillllH]. 

, Although, current humanoid robots more and more resemble human motion, we still need to 

QQ ' emphasize that human joints are (and will probably always remain) significantly more flexible than 

If^ , humanoid robot joints. Namely, each humanoid joint consists of a pair of coupled segments with 

CO ' only Eulerian rotational degrees of freedom. On the other hand, in each human synovial joint, 

ly-^ , besides gross Eulerian rotational movements (roll, pitch and yaw), we also have some hidden and 

' restricted translations along (X,Y, Z)—a.xes. For example, in the knee joint, patella (knee cap) 

moves for about 7-10 cm from maximal extension to maximal flexion). It is well-known that even 
greater are translational amplitudes in the shoulder joint. In other words, within the realm of rigid 
body mechanics, a segment of a human arm or leg is not properly represented as a rigid body fixed 
at a certain point, but rather as a rigid body hanging on rope-like ligaments. More generally, the 
whole skeleton mechanically represents a system of flexibly coupled rigid bodies, technically an 
' anthropomorphic topological product of SE(3)-groups. This implies the more complex kinematics, 

dynamics and control than in the case of humanoid robots i9j . 

This paper presents the scientific body of knowledge behind the sophisticated human motion 
simulator, formulated in the fashion of nonlinear humanoid robotics, called the Human Biodynam- 
ics Engine (HBE), designed over the last five years by the present author at Defence Science & 
Technology Organisation, Australia. The HBE is a sophisticated human neuro-musculo-skeletal dy- 
namics simulator, based on generalized Lagrangian and Hamiltonian mechanics and Lie-derivative 
nonlinear control. It includes 270 active degrees of freedom (DOF), without fingers: 135 rota- 
tional DOF are considered active, and 135 translational DOF are considered passive. The HBE 
incorporates both forward and inverse dynamics, as well as two neural-like control levels. Ac- 
tive rotational joint dynamics is driven by 270 nonlinear muscular actuators, each with its own 
excitation-contraction dynamics (following traditional Hill-Hatze biomechanical models) . Passive 
translational joint dynamics models visco-elastic properties of inter-vertebral discs, joint tendons 
and muscular ligaments as a nonlinear spring-damper system. The lower neural control level re- 
sembles spinal-reflex positive and negative force feedbacks, resembling stretch and Golgi reflexes, 
respectively. The higher neural control level mimics cerebellum postural stabilization and velocity 
target-tracking control. The HBE's core is the full spine simulator, considering human spine as 
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a chain of 26 flexibly-coupled rigid bodies (formally, the product of 26 SE(3)-groups). The HBE 
includes over 3000 body parameters, all derived from individual user data, using standard biome- 
chanical tables. The HBE incorporates a new theory of soft neuro-musculo-skeletal injuries, based 
on the concept of the local rotational and translational jolts, which are the time rates of change of 
the total forces and torques localized in each joint at a particular time instant. 



2 Geometrical Formalism of Human— Robot Biodynamics 
2.1 Configuration Manifold of Idealistic Robot Motion 

Representation of an ideal humanoid-robot motion is rigorously defined in terms of rotational 
constrained S'0(3)-groups [6j[8l[22] in all main robot joints (see Figure [T]). Therefore, the config- 
uration manifold Qrob for humanoid dynamics is defined as a topological product of all included 
S0{?>) groups, Qrob — Y\iSO{'iy . Consequently, the natural stage for autonomous Lagrangian 
dynamics of robot motion is the tangent bundle TQrofuL IS] , and for the corresponding autonomous 
Hamiltonian dynamics is the cotangent bundle r*Qro(D [11 [3]. 
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Figure 1 : The configuration manifold Qrob of the humanoid-robot body is defined as a topological 
product of constrained 50(3) groups, Qrob — Hi SO{iy. 



More precisely, the three-axial 50(3)— group of humanoid-robot joint rotations depends on 
three parameters, Euler joint angles = {ip,ip,9), defining the rotations about the Cartesian 

^In mechanics, to each n— dimensional (nD) configuration manifold Q there is associated its 2nD velocity phase- 
space manifold, denoted by TM and called the tangent bundle of Q. The original smooth manifold Q is called the 
base of TM. There is an onto map n : TM — > Q, called the projection. Above each point x a Q there is a tangent 
space TxQ = ■n~^{x) to Q at x, which is called a fibre. The fibre TxQ C TM is the subset of TM , such that the 
total tangent bundle, TM = | | TxQ, is a disjoint union of tangent spaces TxQ to Q for all points x & Q. From 

mSQ 

dynamical perspective, the most important quantity in the tangent bundle concept is the smooth map d ; Q — > TM , 
which is an inverse to the projection tt, i.e, tt o d = Mq, 'jt{v{x)) = x. It is called the velocity vector-field. Its 
graph {x, v{x)) represents the cross-section of the tangent bundle TM . This explains the dynamical term velocity 
phase-space, given to the tangent bundle TM of the manifold Q. The tangent bundle is where tangent vectors live, 
and is itself a smooth manifold. Vector-fields are cross-sections of the tangent bundle. 

System's Lagrangian (energy function) is a natural energy function on the tangent bundle. 

dual notion to the tangent space TmQ to a smooth manifold Q at a point m is its cotangent space T^Q at 
the same point m. Similarly to the tangent bundle, for a smooth manifold Q of dimension n, its cotangent bundle 
T*Q is the disjoint union of all its cotangent spaces T^Q at all points m a Q, i.e., T*Q = | | T^Q. Therefore, 

the cotangent bundle of an n— manifold Q is the vector bundle T*Q = (TM)* , the (real) dual of the tangent bundle 
TM. The cotangent bundle is where 1-forms live, and is itself a smooth manifold. Covector— fields (1-forms) are 
cross-sections of the cotangent bundle. The Hamiltonian is a natural energy function on the tangent bundle. 
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coordinate triedar (x,y,z) placed at the joint pivot point. Each of the Euler angles are defined in 
the constrained range (— tt, tt), so the joint group space is a constrained sphere of radius tt 6, 8, 2^]. 

Let G = SO{3) = {Ae TWsxsW : A* A = /3,det(A) = 1} be the group of rotations in jt 
is a Lie group and dim(G') — 3. Let us isolate its one-parameter joint subgroups, i.e., consider the 
three operators of the finite joint rotations R^p, R^, Re G 5*0(3), given by 
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corresponding respectively to rotations about x— axis by an angle Lp, about j/— axis by an angle tp, 
and about z— axis by an angle 9. 

The total three-axial joint rotation A is defined as the product of above one-parameter rotations 
R^, RjI,, Rg, i.e., A — R^- R^ ■ Re is equaH 
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cos lp cos ip — COS 9 sin ip sin -0 cos ip cos ip + cos 9 cos p sin ip sin 9 sin ip 
— sin lp cos iy9 — cos 9 sin 1^9 sin ip — sin sin ip + cos cos ip cos ■0 sin cos ip 
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However, the order of these matrix products matters: different order products give different results, 
as the matrix product is noncommutative product. This is the reason why Hamilton's quaternionqj 
are today commonly used to parameterize the 50(3)— group, especially in the field of 3D computer 
graphics. 

/ 1 \ 

The one-parameter rotations i?^, i?^, i?e define curves in 50(3) starting from /3 = 10 

V 1 / 

Their derivatives in = O,'0 = and 9 — belong to the associated tangent Lie algebra so(3). 
That is the corresponding infinitesimal generators of joint rotations - joint angular velocities 
Vip, v^,ve G 50(3) - are respectively given by 
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Moreover, the elements are linearly independent and so 
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The Lie algebra so (3) is identified with 

h 

V G SO (3) given by v 





^ by associating to each v — (u^ , v^,ve) £ 
Then we have the following identities: 



the matrix 



1. u y. V = [u,v\\ and 

2. u ■ V = —\ Tr(?i • v). 



^Note that this product is noncommutative, so it really depends on the order of multiplications. 
■^Recall that the set of Hamilton's quaternions H represents an extension of the set of complex numbers C. We 
can compute a rotation about the unit vector, u by an angle d. The quaternion q that computes this rotation is 



cos — , u sm — 
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The exponential map exp : so (3) 



5*0(3) is given by Rodrigues relation 




, , ^ sni w 
exp(i;) = / + V 

\M\ 

where the norm is given by 

||w|| = ^iv^)^ + iv^)^ + iv3y. 

The dual, cotangent Lie algebra so (3)*, includes the three joint angular momenta Pip,p^,pg S 
so(3)*, derived from the joint velocities v by multiplying them with corresponding moments of 
inertia. 



2.2 Configuration Manifold of Realistic Human Motion 

On the other hand, human joints are more flexible than robot joints. Namely, every rotation in all 
synovial human joints is followed by the corresponding micro-translation, which occurs after the 
rotational amplitude is reached [9]. So, representation of human motion is rigorously defined in 
terms of Euclidean S'i?(3)-groups of full rigid-body motion [TOl [U [H [22] in all main human joints 
(see Figure [2]). Therefore, the configuration manifold Qhum for human dynamics is defined as a 
topological product of all included constrained SE{3) groups, Qrob = YliSE{3y. Consequently, 
the natural stage for autonomous Lagrangian dynamics of human motion is the tangent bundle 
TQhum [5], and for the corresponding autonomous Hamiltonian dynamics is the cotangent bundle 

T*Qw Eiiaiii. 




Figure 2: The configuration manifold Qimm of the human body is defined as a topological product 
of constrained S'i?(3) groups acting in all major (synovial) human joints, Qhum = Ilj SE{3y . 

Briefly, the Euclidean SE(3)-group is defined as a semidirect (noncommutative) product of 
3D rotations and 3D translations, SE{?>) := 50(3) l> R^. Its most important subgroups are the 
following (for technical details see Appendix, as well as p, [TTl [^ ): 
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Subgroup 




Definition 




Sf^(^] PTonD of rotj^tinym 

in 3D (a spherical joint) 


Set of all proper orthogonal 
3 X 3 — rotational matrices 


iS'-£/(2), special Euclidean group 
in 2D (all planar motions) 


Se 


t of all 3 X 3 — matric 
cos 6 sin 9 
— sin 9 cos 9 Vy 
1 


es: 


50(2), group of rotations in 2D 
subgroup of 5'i?(2)-group 
(a revolute joint) 


Set of all proper orthogonal 
2 X 2 — rotational matrices 
included in SE{2) — group 


M"^, group of translations in 3D 
(all spatial displacements) 


Euclidean 3D vector space 



2.3 The Covariant Force Law and Mechanics of Musculoskeletal Injury 

The SE(3)-dynamics applied to human body gives the fundamental law of biomechanics, which is 
the covariant force law [6l [71 [H [22j . It states: 

Force co-vector field = Mass distribution x Acceleration vector field, 

which is formally written (using Einstein's summation convention over repeating indices, with in- 
dices labelling the three Cartesian (X-Y-Z)-translations and the corresponding three Euler angles): 



F„ 



= 1,...,6) 



where denotes the 6 covariant components of the external "pushing" SE(3)-force co- vector field, 
'n^^^v represents the 6x6 covariant components of proximal segment's inertia-metric tensor, while a'^ 
corresponds to the 6 contravariant components of the segment's internal SE(3)-acceleration vector- 
field. This law states that contrary to common perception, acceleration and force are not quantities 
of the same nature: while acceleration is a non-inertial vector field, force is an inertial co-vector 
field. This apparently insignificant difference becomes crucial in injury prediction/prevention, as 
formalized below. Geometrical elaboration of the covariant force law (briefly shown in Figure [3]) 
is fully elaborated in [71 [22] 




Configuration Manifold 



if = if + T'j^v^v'' = x' + r^-j.i-'i 

Levi-Civita Connections 



Figure 3: Riemannian-symplectic geometry of the covariant force law. 

Now we come to injury prediction. It was shown in [51 [7] that the general cause spinal and other 
musculoskeletal injuries is the SE(3)-jolt, which is a sharp and sudden change in the SE(3)-force 
acting on the mass-inertia distribution of the proximal segment to the injured joint. The SE(3)-jolt 
is a 'delta'-change in a total 3D force-vector acting on joint coupled to a total 3D torque- vector. 
In other words, the SE(3)-jolt is a sudden, sharp and discontinues shock in all 6 coupled DOF, 
distributed along the three Cartesian {x, y, z)-translations and the three corresponding Euler angles 
around the Cartesian axes: roll, pitch and yaw. The SE(3)-jolt is rigorously defined in terms of 
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differential geometry [Sld^]- Briefly, it is the absolute time-derivative of the covariant force 1-form 
acting on the joint. 

Formally, the covariant (absolute, Bianchi) time-derivative ^(•) of the covariant SE(3)-force 
defines the corresponding external "striking" SE(3)-jolt co-vector field: 



where ^{a'^) denotes the 6 contravariant components of the proximal segment's internal SE(3)-jerk 
vector- field and overdot (') denotes the time derivative. Fj^^ are the ChristofFel's symbols of the 
Levi-Civita connection for the SE(3)-group, which are zero in case of pure Cartesian translations 
and nonzero in case of rotations as well as in the full-coupling of translations and rotations. 

In particular, the spine, or vertebral column, dynamically represents a chain of 26 movable ver- 
tebral bodies, joint together by transversal viscoelastic intervertebral discs and longitudinal elastic 
tendons (see Figure [3]). Textbooks on functional anatomy describe the following spinal movements: 
(a) cervical intervertebral joints allow all three types of movements: flexion and extension (in the 
sagittal plane), lateral flexion (in the frontal plane) and rotation (in the transverse plane); (b) tho- 
racic joints allow rotation and lateral flexion (limited by ribs), while flexion/extension is prevented; 
and (c) lumbar joints allow flexion/extension as well as limited lateral flexion, while rotation is 
prevented. This popular picture is fine for the description of safe spinal movements; however, to be 
able to predict and prevent spinal injuries (both soft ones related to the back-pain syndrome and 
hard ones related to discus hernia), which are in the domain of unsafe intervertebral movements, 
a much more rigorous description is needed. The main cause of spinal injuries is the SE(3)-jolt, a 
shock that breaks the spinal structure and/or function. 




Figure 4: Reference frame of the Human Biodynamics Engine (HBE). The purpose of the HBE simulator 
is prediction of the risk of soft spinal and other musculo-skeletal injuries, as well as biodynamical behavior 
modelling. 

2.4 Lagrangian Formulation of Biodynamics 

The general form of Lagrangian human/humanoid biodynamics on the corresponding Riemannian 
tangent bundles TQrob and TQhum of the configuration manifolds Qrob and Qhum (precisely derived 
in OIHIH]) can be formulated in a unified form as: 



where n denotes the number of DOF for both rihum and rirob, L = L{t, x, x) : TQ — > R is 
the human/humanoid Lagrangian function, defined on the {2n + l)-dimensional jet manifolds^ 

^In mechanics, we consider a pair of maps /i, /2 : M — > Q from the real line R, representing the time t— axis, into 
a smooth nD configuration manifold Q. We say that the two maps fi = fi{t) and /2 = /2(f) have the same fc— jet 




(1) 



7, -^x^ — (^7 -^T •^) : 



{i = 1, ...,n) 



(2) 
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Xrob = JrobO^^ Qrob) = M X TQrob and Xhuni = Jkunii^^ Qhuva) = M X TQhuni, respectively, with 
local canonical variables {t; a;*^^; iro&) ^hum'^^hum)^ respectively. Its coordinate and velocity 

partial derivatives are respectively denoted by L^i and L^i . 

2.5 Local Muscular Mechanics 

The right-hand side terms J-i{t, x, x) of ([2]) denote any type of external torques and forces, including 
excitation and contraction dynamics of muscular-actuators and rotational dynamics of hybrid 
robot actuators, as well as (nonlinear) dissipative joint torques and forces and external stochastic 
perturbation torques and forces. In particular, we have [SI [Gj El [8] ) : 

1. Synovial joint dynamics, giving the first stabilizing effect to the conservative skeleton dy- 
namics, is described by the (a;,i)-form of the Rayleigh - Van der Pol's dissipation function 

n 

1=1 

where ai and /3j denote dissipation parameters. Its partial derivatives give rise to the viscous- 
damping torques and forces in the joints 

which are linear in x^ and quadratic in a;*. 

2. Muscular dynamics, giving the driving torques and forces jT™"^'^''^ — x, i) with 
(i = 1, . . . ,n) for RHB, describes the internal excitation and contraction dynamics of equivalent 
muscular actuators [T?|. 

(a) Excitation dynamics can be described by an impulse force-time relation 

Fl'^P = - e-*/"') if stimulation >0 

Fl'^P = F°e~'/^' if stimulation =0, 

where denote the maximal isometric muscular torques and forces, while denote the associated 
time characteristics of particular muscular actuators. This relation represents a solution of the 
Wilkie's muscular active-state element equation [T3] 

/i + 7Ai = 75A, m(0) = 0, 0<S'<1, 

where /i — fj,{t) represents the active state of the muscle, 7 denotes the element gain, A corresponds 
to the maximum tension the element can develop, and S — S{r) is the 'desired' active state as a 
function of the motor unit stimulus rate r. This is the basis for the RHB force controller. 

(b) Contraction dynamics has classically been described by the Hill's hyperbolic force-velocity 
relation [T3] 

pHUi ^ - %a.,xJ) 

{Siji^ + bi) 

where and bi denote the Hill's parameters, corresponding to the energy dissipated during the 
contraction and the phosphagenic energy conversion rate, respectively, while Sij is the Kronecker's 
(5— tensor. 

In this way, RHB describes the excitation/contraction dynamics for the ith equivalent muscle- 
joint actuator, using the simple impulse-hyperbolic product relation 

T^"-'^''{t,x,x) = Fl'^P X Fl^'^K 

f at a specified time instant to € K, iff: 

1. /i(t) = f2{t) at to 6 M; and also 

2. the first k terms of their Taylor-series expansions around to G R are equal. 

The set of all fc-jets f : M Q is the fc-jet manifold J*(R, Q). In particular, ji(R, Q) S R X TQ (for technical 
details, see pi22l'). 
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Now, for the purpose of biomedical engineering and rehabilitation, RHB has developed the 
so-called hybrid rotational actuator. It includes, along with muscular and viscous forces, the D.C. 
motor drives, as used in robotics [TSl [51 [B] 

with 

lkik{t) + Rkik{t) + CkXk{t) = Uk{t), 

where k ~ 1, . . . , n, ik{t) and Uk{t) denote currents and voltages in the rotors of the drives, Ik 
and Ck are resistances, inductances and capacitances in the rotors, respectively, while Jk and Bk 
correspond to inertia moments and viscous dampings of the drives, respectively. 

Finally, to make the model more realistic, we need to add some stochastic torques and forces 

mm 

where Bij [x{t), t] represents continuous stochastic diffusion fluctuations, and (t) is an A^— variable 
Wiener process (i.e. generalized Brownian motion), with dW^{t) — W^{t + dt) — W^{t) for 
j = l,...,N. 

2.6 Hamiltonian Biodynamics and Reflex Servo— Control 

The general form of Hamiltonian human/humanoid biodynamics on the corresponding symplectic 
cotangent bundles T*Qrob and T*Qhum of the configuration manifolds Qrob and Qhum (derived in 
imSllS]) is based on the affine Hamiltonian function Ha : T*Q R, in local canonical coordinates 
on T*Q given as 

Ha(x,p,u) = Ho{x,p) - W{x,p)uj, (3) 

where Ho{x,p) is the physical Hamiltonian (kinetic + potential energy) dependent on joint co- 
ordinates x^ and canonical momenta p\ = H^(x,p), [j = 1, . . . , m < n are the coupling 
Hamiltonians corresponding to the system's active joints and Ui = Ui{t, x,p) are (reflex) feedback- 
controls. Using ([3]) we come to the affine Hamiltonian control HBE-system, in deterministic form 

= dp^ Ha - W uj +dp^R, (4) 
Pi =Ti- d^iHa + d^iW Uj + d^iR, 
= -du^Ha = H\ 

X\0)^XI,, P^iO)=pl 

{i^l,...,n; j = 1,. . . , Q < n), 

(where 9„ = d/du, Ti = Ti{t,x,p), Hq = Ho{x,p), W = H^{x,p), Ha = Ha{x,p, u), R = R{x,p)), 
as well as in the fuzzy-stochastic form [TJ [23] 

dq' = (dp^Hoicr^) - dp^W{ap) uj + dp^R) dt, 

dp, = B,,[x\t),t]dW^t) + (5) 

{Ti - d^iHoiafj,) + d^,H^{af,) uj + d^iR) dt, 
do' = -daMa^) dt = W{a^) dt, 
xX0)=4, K(0)=p° 

In ([I])-©, R = R{x,p) denotes the joint (nonlinear) dissipation function, o* are affine system 
outputs (which can be different from joint coordinates); {cr};^ (with /i > 1) denote fuzzy sets 
of conservative parameters (segment lengths, masses and moments of inertia), dissipative joint 
dampings and actuator parameters (amplitudes and frequencies), while the bar (.) over a variable 
denotes the corresponding fuzzified variable; Bij\q^(t),t] denote diffusion fluctuations and W^{t) 
are discontinuous jumps as the n-dimensional Wiener process. 

In this way, the force HBE servo-controller is formulated as affine control Hamiltonian-systems 
(jlHS]), which resemble an autogenetic motor servo [16], acting on the spinal-reflex level of the 
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human locomotion control. A voluntary contraction force F of human skeletal muscle is reflexly 
excited (positive feedback +F^^) by the responses of its spindle receptors to stretch and is reflexly 
inhibited (negative feedback ~F^^) by the responses of its Golgi tendon organs to contraction. 
Stretch and unloading reflexes are mediated by combined actions of several autogenetic neural 
pathways, forming the so-called 'motor servo.' The term 'autogenetic' means that the stimulus 
excites receptors located in the same muscle that is the target of the reflex response. The most 
important of these muscle receptors are the primary and secondary endings in the muscle-spindles, 
which are sensitive to length change - positive length feedback and the Golgi tendon organs, 

which are sensitive to contractile force - negative force feedback —F^^. 

The gain G of the length feedback +F~^ can be expressed as the positional stiffness (the ratio 
G oi S — dF/dx of the force-F change to the length-x change) of the muscle system. The greater 
the stiffness 5', the less the muscle will be disturbed by a change in load. The autogenetic circuits 
-\-F~^ and —F~^ appear to function as servoregulatory loops that convey continuously graded 
amounts of excitation and inhibition to the large (alpha) skeletomotor neurons. Small (gamma) 
fusimotor neurons innervate the contractile poles of muscle spindles and function to modulate 
spindle-receptor discharge. 



2.7 Cerebellum-Like Velocity and Jerk Control 

Nonlinear velocity and jerk (time derivative of acceleration) servo-controllers [21], developed in 
HBE using the Lie-derivative formalism, resemble self-stabilizing and adaptive tracking action of 
the cerebellum jl7 |. By introducing the vector-fields / and g, given respectively by 

/ = (a^^ilo, -a,,i7o) , g = {-dpM\ d,.W) , 

we obtain the affine controller in the standard nonlinear MIMO-system form (see [l8l fT9l [8]) 

x^^ ,f{x)+g{x)uj. (6) 



Finally, using the Lie derivative formalism [22l 125)^1 and applying the constant relative degree r 
to all HE joints, the control law for asymptotic tracking of the reference outputs o-^ = could 
be formulated as (generalized from [TB]) 



(7) 



where c^-i are the coefficients of the linear differential equation of order r for the error function 
e{t) = x^{t)-o'„{t) 

e('') + c^„ie('^-i) + • • • + cie(i) + coe = 0. 

The affine nonlinear MIMO control system ^ with the Lie-derivative control law (O resembles 
the self-stabilizing and synergistic output tracking action of the human cerebellum [531 dS] . To 



^Let F(M) denote the set of all smooth (i.e., C°°) real valued functions / : M — » K on a smooth manifold M, 
V{M) — the set of all smooth vector-fields on M, and V*{M) - the set of all differential one-forms on M . Also, let 
the vector— field C, G V{M) be given with its local flow (p^. : M ^ M such that at a point x G M, ^|t=o <j>t^ = CC^^), 
and 0* representing the pull-back by <f>f The Lie derivative differential operator £^ is defined: 

(i) on a function / g F{M) as 

£ :F(M)^F(M), Ccf = ^Mf)\t=o, 

at 

(ii) on a vector-field r] G V{M) as 

£ : viM) ^ ViM), Cell = ^U*tri)\t=o = [Cv] 

at 

- the Lie bracket, and 

(iii) on a one-form a G y*(M) as 

£^:l/*(M)^y*(M), = ^('^'t")lt=0- 

In general, for any smooth tensor field T on M, the Lie derivative geometrically represents a directional 

derivative of T along the flow (p^. 
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make it adaptive (and thus more realistic), instead of the 'rigid' controUer ([7]), we can use the 
adaptive Lie-derivative controller, as explained in the seminal paper on geometrical nonlinear 
control gO]. 

2.8 Cortical— Like Fuzzy— Topological Control 

For the purpose of our cortical control, the dominant, rotational part of the human configuration 
manifold M^, could be first, reduced to an iV-torus, and second, transformed to an A''-cube 
('hyper-joystick'), using the following topological techniques (see |5[ 1^ [^ ). 

Let denote the constrained unit circle in the complex plane, which is an Abelian Lie group. 
Firstly, we propose two reduction homeomorphisms, using the Cartesian product of the constrained 
S'0(2)-groups: 

S0{3) w 50(2) X SO{2) x S0{2) and SO{2) « S^. 

Next, let be the unit cube [0, 1]^ in and an equivalence relation on obtained 
by 'gluing' together the opposite sides of , preserving their orientation. Therefore, can be 
represented as the quotient space of M.^ by the space of the integral lattice points in M^, that is 
an oriented and constrained A^-dimensional torus : 

N 

R^/Z^ ~ n -^^^ = « = 1, • ■ • , ^) : mod27r} = . (8) 

i=l 

Its Euler-Poincare characteristic is (by the De Rham theorem) both for the configuration manifold 
and its momentum phase-space T*T^ given by (see [22) ) 

N 

where bj, are the Betti numbers defined as 

bP^ (^^^,...6^-1 = TV, 
iO<p<N). 

Conversely by 'ungluing' the configuration space we obtain the primary unit cube. Let '^*^ 
denote an equivalent decomposition or 'ungluing' relation. According to Tychonoff's product- 
topology theorem [8l[22], for every such quotient space there exists a 'selector' such that their 
quotient models are homeomorphic, that is, / / ~*. Therefore represents a 'selector' 

for the configuration torus and can be used as an iV-directional '^-command-space' for the 
feedback control (FC). Any subset of degrees of freedom on the configuration torus representing 
the joints included in HB has its simple, rectangular image in the rectified q-command space - 
selector , and any joint angle has its rectified image q*. 

In the case of an end-effector, reduces to the position vector in external-Cartesian coordinates 
z"" (r = 1, . . . , 3). If orientation of the end-effector can be neglected, this gives a topological solution 
to the standard inverse kinematics problem. 

Analogously, all momenta pi have their images as rectified momenta pi in the j5-command space 
- selector 1^ . Therefore, the total momentum phase-space manifold T*T^ obtains its 'cortical 

image' as the (q,p)-command space, a trivial 2iV-dimensional bundle x . 

Now, the simplest way to perform the feedback FC on the cortical (q,p)-command space 
X , and also to mimic the cortical-like behavior, is to use the 2N~ dimensional fuzzy-logic 

controller, in much the same way as in the popular 'inverted pendulum' examples (see |26j). 

We propose the fuzzy feedback-control map S that maps all the rectified joint angles and 

momenta into the feedback-control one-forms 

R:{q'{t),p,{t))^u,{t,q,p), (9) 



b^ =N,... 

b^ = 1, 
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so that their correspondmg universes of discourse, — (Qmax ~ 'imin)^ — {p''^^^ — p™™) and 
Ui = {u^"-^ — M™™), respectively, are mapped as 

AT AT TV 

^ WQ' ^\{P^^X{lJ^■ (10) 

i—1 i—1 i—1 

The 27V-dimensional map S (|9ll0p represents a fuzzy inference system, defined by (adapted 
from 

1. Fuzzification of the crisp rectified and discretized angles, momenta and controls using Gaussian- 
bell membership functions 

/.,(X) = (fc = l,2,...,9), 

Z<7k 

where x G Z? is the common symbol for g*, pi and Ui{q,p) and D is the common symbol for 
Q\ Pi and i; the mean values nik of the nine partitions of each universe of discourse D are 
defined as = XkD+Xmim with partition coefficients Xk uniformly spanning the range of D, 
corresponding to the set of nine linguistic variables L = {NL, NB, NAI, NS, ZE, PS, PM, 
PB,PL}: standard deviations are kept constant ak = D/9. Using the linguistic vector L, 
the 9x9 FAM (fuzzy associative memory) matrix (a 'linguistic phase-plane'), is hcuristically 
defined for each human joint, in a symmetrical weighted form 

Mfei = ■^ki exp{-5Q[\k + u{q,p)]'^}, (fc, I = 1, 9) 

with weights ruki S {0.6, 0.6, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 1.0}. 

2. Mamdani inference is used on each FAM-matrix /x^,; for all human joints: 

(i) fJ,{q^) and /i(pi) are combined inside the fuzzy IF-THEN rules using AND (Intersection, 
or Minimum) operator, 

(^kiuiiq,p)] = min{/i;,;(g*), ^ikliPi)}■ 

(ii) the output sets from different IF-THEN rules are then combined using OR (Union, or 
Maximum) operator, to get the final output, fuzzy-covariant torques, 

fj,[u,{q,p)] ^m&x{p^[u,{q,p)]}. 

k 

3. Defuzzification of the fuzzy controls iJ,[ui{q,p)] with the 'center of gravity' method 

/ ^ J f^[ui{q,p)]dui 

Ui[q,p) = TT- , 

J dui 

to update the crisp feedback-control one-forms Ui = Ui{t, q,p). 

Now, it is easy to make this top-level controller adaptive, simply by weighting both the above 
fuzzy-rules and membership functions, by the use of any standard competitive neural-network 
(see, e.g., [H]). Operationally, the construction of the cortical (g,p)-comniand space x 1^ and 
the 2iV-dimensional feedback map S (|9ll0p , mimic the regulation of the motor conditioned reflexes 
by the motor cortex fVP . 



3 HBE Simulation Examples 

The first version of the HBE simulator had the full human-like skeleton, driven by the generalized 
Hamiltonian dynamics (including muscular force-velocity and force-time curves) and two levels 
of reflex-like motor control (simulated using the Lie derivative formalism) [31 [51 [7]. It had 135 
purely rotational DOF, strictly following Figure [H It was created for prediction and prevention of 
musculo-skeletal injuries occurring in the joints, mostly spinal (intervertebral, see Figure [3]). Its 
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Figure 5: Sample output from the Human Biodynamics Engine: running simulation with the speed of 6 
m/s - 3D animation view-port. 



performance looked kinematically realistic, while it was not possible to validate the driving torques. 
It included a small library of target movements which were followed by the HBE's Lie-derivative 
controllers with efficiency of about 90% (see Figures |7| and [5]) . 

The HBE also includes a generic crash simulator, based on the simplified road- vehicle impact 
simulator (see Figure [9]). While implementing the generic crash simulator, it became clear that 
purely rotational joint dynamics would not be sufficient for the realistic prediction of musculo- 
skeletal injuries. Especially, to simulate the action of a Russian aircraft ejection-seat currently 
used by the American Space-shuttle, we needed, strictly following Figure |2l to implement micro- 
translations in the intervertebral joints (see Figures fTUl and [TT|) . as the seat provides the full body 
restraint and the ejection rockets, firing with 15 g for 0.15 s - can only compress the spine, without 
any bending at all. 

In this way a full rotational -I- translational biodynamics simulator has been created with 270 
DOF in total (not representing separate fingers). The 'HBE-simulator' has been kinematically 
validated [55] against the standard biomechanical gait-analysis system 'Vicon' 

4 Complexity of Biodynamical Behavior 
4.1 Biodynamical 'Entanglement' 

From the standard engineering viewpoint, having two systems combined, in the case of biodynamics 
- biological and mechanical, as a single 'working machine', we can expect that the total 'machine' 
complexity equals the sum of the two partial ones. For example, electrical circuitry has been a 
standard modelling framework in neurophysiology (A. Hodkgin and A. Huxley won a Nobel Prize 
for their circuit model of a single neuron, the celebrated HH-neuron ,34j)- Using the HH-approach 
for modelling human neuro-muscular circuitry as electrical circuitry, we get an electro-mechanical 
model for our bio-mechanical system, in which the superposition of complexities is clearly valid. 

On the other hand, in a recent research on dissipative quantum brain modelling, one of the 
most popular issues has been quantum entanglemeniu between the brain and its environment (see 

^Entanglement is a term used in quantum theory to describe the way that particles of energy/matter can become 
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Figure 6: Matching tiie 'Vicon' output witli the 'HBE' output for the riglit-hip angular velocity around 
the dominant X-axis, while walking on a treadmill with a speed of 4 km/h. 

[30l[3T]) where the brain-environment system has an entangled 'memory' state (identified with its 
ground state), that cannot be factorized into two single-mode statesjfl Similar to this microscopic 
brain-environment entanglement, we conjecture the existence of a macroscopic neuro-mechanical 
entanglement between the operating modes of our neuro-muscular controller and purely mechanical 
skeleton (see [35]). 

In other words, we suggest that the diffeomorphism between the brain motion manifold (iV— cube) 
and the body motion manifold (which can be reduced to the constrained TV— torus), described 
as the cortical motion control, can be considered a 'long-range correlation', thus manifesting the 
'biodynamical entanglement'. 

4.2 Biodynamical Self— Assembly 

In the framework of human motion dynamics, self-assembly represents adaptive motor control, i.e., 
physiological motor training performed by iteration of conditioned reflexes. For this, a combination 
of supervised and reinforcement training is commonly used, in which a number of (nonlinear) 
control parameters are iteratively adjusted similar to the weights in neural networks, using either 
backpropagation-type or Hebbian-type learning, respectively (see, e.g., [IS]). Every human motor 
skill is mastered using this general method. Once it is mastered, it becomes smooth and energy- 
efficient, in accordance with Bernstein's motor coordination and dexterity (see [3HI31]). 

Therefore, biodynamical self-assembly clearly represents an 'evolution' in the parameter-space 
of human motion control. One might argue that such an evolution can be modelled using CA. 
However, this parameter-space, though being a dynamical and possibly even a contractible struc- 
ture, is not an independent set of parameters - it is necessarily coupled to the mechanical skeleton 
configuration space, the plant to be controlled. 

The system of 200 bones and 600 muscles can an produce infinite number of different move- 
ments. In other words, the output-space dimension of a skilled human motion dynamics equals 
infinity - there is no upper limit to the number of possible different human movements, starting 
with simple walk, run, jump, throw, play, etc. Even for the simplest motions, like walking, a child 
needs about 12 months to master it (and Honda robots took a decade to achieve this). 

Furthermore, as human motion represents a simplest and yet well-defined example of a general 
human behavior, it is possible that other human behavioral and performance skills are mastered 
(i.e., self-assembled) in a similar way. 

correlated to predictably interact with each other regardless of how far apart they are; this is called a 'long-range 
correlation'. 

*In the Vitiello— Pessa dissipative quantum brain model I30II31I . the evolution of a memory system was represented 
as a trajectory of given initial condition running over time-dependent states, each one minimizing the free energy 
functional. 
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Figure 7: The HBE simulating a jump-kick: a 3D viewer. 



4.3 Biodynamical Synchronization 

The route to simpUcity in biodynamics is synchronization. Both synchronization and phase-locking 
are ubiquitous in nature as well as in human brain (see [35[ I36[ I37[ I38j ) . Synchronization can occur 
in cyclic forms of human motion (e.g., walking, running, cycling, swimming), both externally, in 
the form of oscillatory dynamics, and internally, in the form of oscillatory cortical-control. This 
oscillatory synchronization, e.g., in walking dynamics, has three possible forms: in-phase, anti- 
phase, and out-of-phase. The underlying phase-locking properties determined by type of oscillator 
(e.g., periodic/chaotic, relaxation, bursting, pulse-coupled, slowly connected, or connections with 
time delay) involved in the cortical control system (motion planner). According to Izhikevich- 
Hoppensteadt work (ibid) , phase-locking is prominent in the brain: it frequently results in coherent 
activity of neurons and neuronal groups, as seen in recordings of local field potentials and EEG. 
In essence, the purpose of brain control of human motion is reduction of mechanical configuration 
space; brain achieves this through synchronization. 

While cyclic movements indeed present a natural route to Biodynamical synchronization, both 

^Periodic bursting behavior in neurons is a recurrent transition between a quiescent state and a state of repetitive 
firing. Three main types of neural bursters are: (i) paraboUc bursting ('circle/circle'), (ii) square-wave bursting 
('fold/homochnic'), and (iii) elliptic bursting ('subHopf/fold cycle'). Most burster models can be written in the 
singularly perturbed form: 

^ = fix,y), y = ^ig{x,y), 

where x G M"" is a vector of fast variables responsible for repetitive firing (e.g., the membrane voltage and fast 
currents). The vector y G is a vector of slow variables that modulates the firing (e.g., slow (in) activation dynamics 
and changes in intracellular Ca^+ concentration). The small parameter ^ << 1 is a ratio of fast/slow time scales. 
The synchronization dynamics between bursters depends crucially on their spiking frequencies, i.e., the interactions 
are most effective when the presynaptic inter-spike frequency matches the frequency of postsynaptic oscillations. 
The synchronization dynamics between bursters in the cortical motion planner induces synchronization dynamics 
between upper and lower limbs in oscillatory motions. 
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Figure 8: The HBE simulating a jump-kick: calculating joint angles and muscular torques. 

on the dynamical and cortical-control level, the various forms of synchronized group behavior in 
sport (such as synchronized swimming, diving, acrobatics) or in military performance represent 
the imperfect products of hard training. The synchronized team performance is achievable, but 
the cost is a difficult long-term training and sacrifice of one's natural characteristics. 

5 Conclusion 

We have presented various aspects of development of the Human Biodynamics Engine. The HBE 
geometry is based on anthropomorphic tree of Euclidean SE(3)-groups. Its dynamics was at 
first Lagrangian and later changed to Hamiltonian (dynamically equivalent, but superior for con- 
trol). Its actuators are 'equivalent muscles', following classical Hill-Hatze muscular mechanics. Its 
reflexes follow Houk's 'autogenetic' stretch-Golgi prescription. Its 'cerebellum' is modelled using 
Lie-derivative formalism. Its brain is fuzzy-topological. Its complexity shows biodynamical 'entan- 
glement', self-assembly and oscillatory synchronization. Its simulations demonstrate the necessity 
of micro-translations in the human joints, which cannot exist in robots. The main purpose for its 
development has been prediction of neuro-musculo-skeletal injuries. For this purpose, the concept 
of rotational (soft) and translational (hard) jolts has been developed and implemented in HBE. 
The HBE simulator is currently under the thorough validation process. Kinematic validation has 
mostly been completed, while for the validation of torques and forces we are still lacking adequate 
in vivo measurement technology. 
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Figure 9: The HBE simulating the frontal road- vehicle crash into the fixed wall with a speed of 70 
km/h: before the impact (up) and 0.12 s after the impact. 



6 Appendix: The S* £"(3) — Group of General Rigid Motions 

The special Euclidean group SE{3) :~ SO{3)l>M.^ , (the semidirect product of the group of rotations 
with the corresponding group of translations), is the Lie group consisting of isometrics of the 
Euclidean 3D space (see PfTT1[22]). 

An element of SE{3) is a pair {A, a) where A e SO{S) and aeR^. The action of SE{3) on 
is the rotation A followed by translation by the vector a and has the expression 

(A, a) ■ X — Ax + a. 

The Lie algebra of the Euclidean group SE{3) is se(3) x ^itj^ tj^g Lie bracket 

[(f , u), {ri, v)]^{Cxr),^xv-'ijxu). (11) 
Using homogeneous coordinates, we can represent SE{3) as follows, 

SE{3)^ ^ ^ eGL(4,M) :i?eS'0(3),p6lE 

with the action on given by the usual matrix-vector product when we identify R^ with the 
section R^ x {1} c M'^. In particular, given 



R p 
1 



e SE{3), 
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Figure 10: The HBE simulating the effect of an aircraft pilot-seat ejection to human spine com- 
pression: before the seat ejection (left) and after ejection (right). 



and g G K."^, we have 

or as a matrix- vector product, 



g-q = Rq+p, 



R p 
1 



Rq+p 
1 



The Lie algebra of SE{3), denoted 5e(3), is given by 



se(3) 



CO V 





e M4(M) : uj e 5o(3), w e mH , 



where the attitude (or, angular velocity) matrix lu 



5o(3) is given by 



-LUz UJy 

—LOy LOx 



The so-called exponential map, exp : se(3) SE{3), is given by 



where 



exp 



A = I 












J 





1 — cos 

and exp(ti;) is given by the Rodriguez' formula, 



, , ^ sm w 
expH = 



exp(a;) Av 
1 



kj - sm LU 



1 — cos 



L0\ 
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Time [s] 

Figure 11: The HBE calculating translational forces distributed along the spinal joints during the 
seat ejection. 
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' Abstract 

' This paper presents the scientific body of knowledge behind the Human Biodynamics En- 

gine (HBE), a human motion simulator developed on the concept of Euclidean motion group 
<—{ • SE(3), with 270 active degrees of freedom, force-velocity-time muscular mechanics and two- 

h—i ' level neural control - formulated in the fashion of nonlinear humanoid robotics. The following 

QQ , aspects of the HBE development are described: geometrical, dynamical, control, physiological, 

' AI, behavioral and complexity, together with several simulation examples. 

I I , Index Terms: Human Biodynamics Engine, Euclidean SE(3)-group, Lagrangian/Hamiltonian 

' biodynamics. Lie-derivative control, muscular mechanics, fuzzy-topological coordination, 

^) \ biodynamical complexity, validation, application 

.2 ■ 

X) ! 1 Introduction 

rV' 

I i ' Both human biodynamics and humanoid robotics are devoted to studying human-like motion. 

They are both governed by Newtonian dynamical laws and reflex-like nonlinear controls [?, ?, ?, 
■ ? ? ? ? ?i 

, Although, current humanoid robots more and more resemble human motion, we still need to 

QQ ' emphasize that human joints are (and will probably always remain) significantly more flexible than 

If^ , humanoid robot joints. Namely, each humanoid joint consists of a pair of coupled segments with 

CO ' only Eulerian rotational degrees of freedom. On the other hand, in each human synovial joint, 

ly-^ , besides gross Eulerian rotational movements (roll, pitch and yaw), we also have some hidden and 

' restricted translations along (X,Y, Z)—axes. For example, in the knee joint, patella (knee cap) 

moves for about 7-10 cm from maximal extension to maximal flexion). It is well-known that even 
greater are translational amplitudes in the shoulder joint. In other words, within the realm of rigid 
body mechanics, a segment of a human arm or leg is not properly represented as a rigid body fixed 
at a certain point, but rather as a rigid body hanging on rope-like ligaments. More generally, the 
whole skeleton mechanically represents a system of flexibly coupled rigid bodies, technically an 
' anthropomorphic topological product of SE(3)-groups. This implies the more complex kinematics, 

dynamics and control than in the case of humanoid robots [?] . 

This paper presents the scientific body of knowledge behind the sophisticated human motion 
simulator, formulated in the fashion of nonlinear humanoid robotics, called the Human Biodynam- 
ics Engine (HBE), designed over the last five years by the present author at Defence Science & 
Technology Organisation, Australia. The HBE is a sophisticated human neuro-musculo-skeletal dy- 
namics simulator, based on generalized Lagrangian and Hamiltonian mechanics and Lie-derivative 
nonlinear control. It includes 270 active degrees of freedom (DOF), without fingers: 135 rota- 
tional DOF are considered active, and 135 translational DOF are considered passive. The HBE 
incorporates both forward and inverse dynamics, as well as two neural-like control levels. Ac- 
tive rotational joint dynamics is driven by 270 nonlinear muscular actuators, each with its own 
excitation-contraction dynamics (following traditional Hill-Hatze biomechanical models) . Passive 
translational joint dynamics models visco-elastic properties of inter-vertebral discs, joint tendons 
and muscular ligaments as a nonlinear spring-damper system. The lower neural control level re- 
sembles spinal-reflex positive and negative force feedbacks, resembling stretch and Golgi reflexes, 
respectively. The higher neural control level mimics cerebellum postural stabilization and velocity 
target-tracking control. The HBE's core is the full spine simulator, considering human spine as 
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